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Abstract. Statistics of the free volume available to individual particles have 
previously been studied for simple and complex fluids, granular matter, amorphous 
solids, and structural glasses. Minkowski tensors provide a set of shape measures that 
are based on strong mathematical theorems and easily computed for polygonal and 
polyhedral bodies such as free volume cells (Voronoi cells). They characterize the 
local structure beyond the two-point correlation function and are suitable to define 
indices < f3^ b < 1 of local anisotropy. Here, we analyze the statistics of Minkowski 
tensors for configurations of simple liquid models, including the ideal gas (Poisson point 
process), the hard disks and hard spheres ensemble, and the Lennard- Jones fluid. We 
show that Minkowski tensors provide a robust characterization of local anisotropy, 
which ranges from (3® ,b ss 0.3 for vapor phases to f3^ ,b —> 1 for ordered solids. We find 
that for fluids, local anisotropy decreases monotonously with increasing free volume 
and randomness of particle positions. Furthermore, the local anisotropy indices /3°' b 
are sensitive to structural transitions in these simple fluids, as has been previously 
shown in granular systems for the transition from loose to jammed bead packs. 
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1. Introduction 

As new experimental techniques become available, more detailed information on the 
nature of the fluid state can be accessed. New scattering experiments reveal the local 
ordering of fluids [1, 2], and confocal microscopy allows for individual particles to 
be tracked [3]. This richness of new data calls for new methods to characterize the 
morphology of equilibrium and nonequilibrium fluids and related systems. 

While fluids of spherical particles usually are globally isotropic, the local 
environments of the particles may be anisotropic (fig. 1). Local anisotropy in jammed 
bead packs has recently been studied in the framework of Minkowski tensors [4]. Local 
anisotropy may, in fact, play a crucial role both for dynamical and equilibrium properties 
of particle systems. At the very least, measures of local anisotropy can be used as a local 
order parameter to discriminate between the various equilibrium and nonequilibrium 
states of an assembly of particles. This article extends the analysis given in [4] to 
important model fluids, namely the ideal gas (Poisson point process), the hard spheres 
and hard disks ensembles, and Lennard- Jones fluids. We demonstrate that the local 
anisotropy yields distinct evidence of the phase behavior of simple fluids; furthermore, 
we show that Minkowski tensors provide a robust means to its characterization. 

As a measure of local geometry in particulate matter, local anisotropy evidently 
competes with existing measures based on nearest-neighbor analysis, including the 
bond-orientational order parameters [5, 6, 7, 8], the isoperimetric ratio [9, 10] and 
related quantities [11], and tensorial measures such as Edwards' configurational 
tensor [12], quadrons [13], and fabric / texture tensors [14]. These measures are, 
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Figure 1. Histograms of local density (x axis) and Voronoi cell anisotropy index {y 
axis) in a Poisson point pattern, a dilute hard disks fluid (77 = 32%), and hard disks 
configurations below and above the phase transition (77 = 61%, 76%). The shape of 
the distribution changes characteristically with the packing density. 
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while undoubtedly widely and successfully applied, reliant on identifying suitable 
neighborhoods. Neighborhoods, by their very nature, change discontinuously, in 
contrast to the shape measures of Voronoi cells used in the present study. The volume 
distributions of Voronoi cells have been studied [15, 16, 17] and the shape of Voronoi cells 
is their next important property. Scalar measures of shape such as the isoperimetric ratio 
contain some signature of anisotropy; however, they do not discern between the generally 
distinct but often correlated concepts of asphericity (meaning deviation form a spherical 
shape) and anisotropy (meaning elongation or directedness) . It is our hypothesis that 
explicit anisotropy measures such as the Minkowski tensors will turn out similarly useful 
for the identification of phase transitions in fluids and other particulate systems, while 
being more generic and robust in their definition than measures based on neighborhoods. 
Explicit anisotropy measures may also be more easily generalized to aspherical particles, 
relevant to liquid crystalline phases [18, 19, 20]. 

The article is organized as follows: In section 2.1, we restate briefly the definition 
of the Minkowski tensors. Section 2.2 discusses how Minkowski tensors may be used to 
characterize anisotropy. In section 2.3, we use Minkowski tensors of Voronoi tessellations 
to characterize the local anisotropy in a point pattern, such as the configuration of 
particles in a fluid. Sections 3 and 4 gives the results of the anisotropy analysis for 
an ideal gas, described by the Poisson point process, and equilibrium systems of hard 
spheres and disks. Section 5 applies the same analysis to point particles interacting with 
a Lennard- Jones potential. Section 6 shows that Minkowski tensors are both robust and 
sensitive measures of local geometry by applying them to the Einstein solid model. 

2. Anisotropy Characterization by Minkowski Tensors 

2.1. Minkowski Tensors 

The description of shape requires the choice of adequate quantities that may be 
associated with the shape in question. Minkowski functionals, also known as 
Quermafiintegrale or intrinsic volumes, are scalar measures of shape defined in convex 
and integral geometry, based on support measures [21, 22]. Equivalently, they may be 
introduced as curvature-weighted integrals, as in this section. 

Given a compact subset of D-dimensional Euclidean space K C 1R D , called a 
body, with a sufficiently smooth boundary, the Minkowski functionals in two and three 
dimensions are defined as 



where G v {r) is a symmetric combination of the principal curvatures Ki of the boundary 
dK. In two dimensions, G\ = \,Gi = k, (there is only one curvature) and in three 
dimensions G\ — \,Gi — \{k>i + ^2)^3 — ^1^2- The Minkowski functionals are, up 
to conventional prefactors, the area, the circumference, and the Euler number (in 2D) 
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Table 1. A basis for the space of rank-two additive, conditionally continuous, motion- 
covariant functionals in terms of Minkowski functionals W v and Minkowski tensors 
W®' b . E is the unit matrix, n = n(r) the local curvature of the boundary contour dK 
at point r, «i(r) and K2(r) the principal curvatures of the boundary surface dK, n{r) 
the outer boundary normal at point r. 



and the volume, the surface area, the integral mean curvature, and the Euler number 
(in 3D) of the body K. The set of Minkowski functionals has a number of properties 
which make them useful for shape analysis of physical systems [23]. In particular, 
by Hadwiger's characterization theorem [21], the set of Minkowski functionals spans 
the space of additive, conditionally continuous and motion invariant functionals of K. 
For such functionals, the Minkowski functionals contain all the relevant morphological 
information about the shape of K. 

A natural generalization of Minkowski functionals are Minkowski tensors, which 
arise by including tensor products of the position vector r and the outer boundary 
normal n in the integrals in eq. (2). In analogy to the scalar case, the Minkowski tensors 
span the space of additive, conditionally continuous and motion- covariant functionals 
[24]. Knowledge of all the Minkowski tensors provides a more complete morphological 
description of the body K than the scalar properties alone, and allows for tensor- valued 
functionals of K to be computed. The set of Minkowski tensors contains, analogous to 
the scalar case, complete morphological information about the shape of K. A basis of 
the space of tensor- valued functionals of rank two is given in table 1. In addition to the 
scalar Minkowski functionals, four Minkowski tensors (six in 3D) are required to span 
the space. 

For the mathematical details we refer the reader elsewhere [25, 26, 27], and instead 

2 

only give a single example, namely the Minkowski tensor W 1 ' of a planar body K which 
is given by 

W?>\K) = \\ dr( niUl niU2 ) (3) 
JdK \ n 2 ni n 2 n 2 J 

where rii are the components of the outer normal vector at the boundary dK of K, and 
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dr is the infinitesimal line element on the boundary OK. W®' 2 {K) can now be written 
in terms of three independent scalar boundary integrals 

Importantly for applications, the Minkowski tensors of a triangulated body (in the planar 
polygon) can be written as a simple sum over all facets and thus calculated in 
O(N) time. For a polygon, the hj(K) in eq. (4) are given by the three sums 

w = lE L(/) W ( 5 ) 
/ 

where the index / labels the edges of the polygon, and are their lengths and 

outer normals respectively. An analogous formula holds in three dimensions, replacing 
edge lengths by facet areas. Explicit formulae for the fast and simple computation of 
the complete set of Minkowski tensors (W 2 '°, W 2 ' , W^' , W?' 2 in 2D and W 2 ' , W 2 ' , 
W 2 2,0 , W 2 ' , W°' 2 , W°' 2 in 3D) can be found in refs. [25, 26]. 

Similar to the tensor of inertia, some of the Minkowski tensors, W%'°, do not depend 
only on the shape of K but also on its location. A suitable origin has to be chosen for 
computation of the Minkowski tensors for a given body. The tensors including only 
normal vectors, W®' 2 , are translation invariant. 



2.2. Anisotropy Characterization of Shapes 

Having assigned a number of tensorial shape descriptors to each body K, we can now 
define an anisotropy index for each body K and each Minkowski tensor W"' b (K) by the 
dimensionless ratio 

where A min (W), X max (W) are the smallest and largest eigenvalues of the D x D matrix 
W; D the spatial dimension. An anisotropy index of j3 — 1 indicates an isotropic body, 
and values ranging from 1 to correspond to increasing degrees of anisotropy. The 
dimensionless anisotropy index is a pure shape measure; it is invariant under isotropic 
scaling of K, i. e. ffi' b (sK) = p? b (K), for all s > 0. 

Note that, in the framework of this analysis, we consider any body isotropic that 
has Minkowski tensors proportional to the unit matrix. Bodies with n-fold rotation 
symmetries through the origin (n > 3) have degenerate eigenvalues due to the additivity 
of the Minkowski tensors; as a consequence, many well-known highly symmetric bodies 
have (3 = 1. This includes, in two dimensions, the sphere and square, and any 
regular polygon; in three dimensions, sphere and cube, and the fee unit cell (truncated 
octahedron) centered at the origin are examples of isotropic bodies. 
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2. 3. Characterization of Point Patterns 

We characterize local anisotropy of point configurations, representing for example the 
centroids of fluid particles, by computing the anisotropy indices (3 of their Voronoi cells. 
Each point in the configuration, termed germ, is assigned a Voronoi cell that 
consists of all points in ~R D closer to than to any other germ g^) ^ g( l ). From 
any point configuration {<7^}, % — 1, . . . , N we thus obtain a set of N bodies {K^}, 
and consequently, N Minkowski tensors {W®' b (K^)}. For the translation invariant 
Minkowski tensors, the choice of origin is irrelevant; for the others, the germ point g® 
is chosen as the origin to compute the tensor W®' b (K^). Voronoi cells are frequently 
used to define the free volume available to an individual particle in a fluid configuration, 
and in glassy and supercooled states, for example [11]. The same concept has been used 
in granular matter for static packings [28] . The Voronoi tessellations of a set of points 
is computed using readily available software [29]. In contrast to the neighborhood 
topology, the shapes of the Voronoi cells are a continuous function of the point pattern. 

For the remainder of the paper, the index % labels the germ point in the point 
configuration. Since we study point configurations consisting of a large number of 
points, all results will be of statistical nature. Averages over the germs in one or more 
configurations are denoted by //( • ); for example, (i((3®' b (K^)) = • J2i (3^' b (K^) 
is the average anisotropy index of Voronoi cells, with the average taken over all the 
germs. In addition, a( ■ ) denotes the standard deviation; for the observable X^ l \ it is 
computed from a 2 (X^) = (N - l)" 1 ■ Ei(^ W - M* (0 )) 2 - For brevity, we identify 

3. Poisson Point Process (Ideal Gas) 

The simplest model of a fluid is the ideal gas; molecules in the gas do not interact, 
and as a consequence, their locations g^ are uncorrelated. Mathematically, ideal 
gas configurations are modeled by a Poisson point process (PPP) [30]. The only 
free parameter in the Poisson point process is the intensity p, which is equal to the 
expectation value of the number of particles per unit volume. This parameter can be 
removed by a rescaling the coordinates in the system. Since the anisotropy indices 
f3®' b , as defined above, are invariant under scaling of the body, the measured degree of 
anisotropy is independent of the intensity p. 

Voronoi tessellations are difficult to treat analytically. Thus, we estimate the 
parameters of the (3®' b distributions numerically. For this purpose, ten realizations of 
Poisson process on a square [0; l) 2 and on a cube [0; l) 3 are generated by first drawing 
the number of germs from a Poisson distribution with mean 16384 and subsequently 
placing the germs in the square or cube; germ point coordinates are drawn independently 
from uniform distributions in on [0; 1). The Voronoi tessellation is computed for each 
configuration with periodic boundary conditions, and Minkowski tensors and functionals 
are computed for the individual Voronoi cells. 
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Table 2. Scalar Minkowski functional of Voronoi cells induced by a Poisson point 
process, p is the intensity of the point process, T the Gamma function, and D the 
spatial dimension. Exact results are taken from [17], numerical values estimated in the 
present study from 160,000 random Voronoi cells. 
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Figure 2. Distribution of the anisotropy index /3°' 2 in configurations of hard disks 
(left) and spheres (right) of various packing fractions. The Poisson limit rj — > is 
indicated by the bold black line. Also shown is the best fit of a Beta distribution 
oc x m (l — x) n to the Poisson limit. 



Table 2 shows results for the average scalar Minkowski functionals of Poisson- 
Voronoi cells. Analytical expressions are known for Wo and W\, and we find numerical 
values close to the exact results [15, 17]. The authors are not aware of an analytical 
expression for the average integral mean curvature //(W2) of Poisson- Voronoi cells; the 
same is true for the standard deviations of the respective distributions. The scalar 
Minkowski functionals (area, perimeter in 2D and volume, surface area, integral mean 
curvature in 3D) are found to be compatible with generalized Gamma distributions in 
agreement with the literature (see fig. 8 in appendix and refs. [31, 17]). 

The mean and standard deviation of the distribution of the anisotropy indices 
of the Poisson- Voronoi cells are shown in table 3. Since the anisotropy indices are 
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Table 3. Mean and standard deviation of the anisotropy indices /3° ,b of the Voronoi 
cells induced by a Poisson point process. The intensity p of the Poisson point process 
does not affect the dimensionless anisotropy indices. 

invariant under scaling, the length scale p does not enter, and the tabulated values 
are fingerprints of a point pattern with uniformly and independently distributed germ 
points. The full distribution of /3°' 2 values is shown in fig. 2 (bold lines). In the 3D 
case, it is approximately a Beta distribution, which is remarkable because the Beta 
distribution is closely related to the Gamma distribution!. The agreement is less good 
in the 2D case. 

4. Equilibrium Hard Spheres and Disks 

The hard disks and hard spheres model is beyond doubt one of the most important 
generalizations of the ideal gas. This model, consisting in spherical particles without 
interaction other than elastic hard repulsion, possesses a fluid-solid phase transition 
[32, 33] and - in the 3D case - vitreous and jammed states [34]. In contrast to 
the ideal gas discussed above, the hard spheres model does have an intrinsic length 
scale by virtue of the particle radius and packing fraction 77. This work focuses on 
monodisperse particles, even though the method is applicable to more general models 
when one replaces the Voronoi tessellation by a more general, for example the Apollonius 
tessellation [17]. 

Two different procedures have been employed to generate realizations of hard disks 
and hard spheres configurations: First, Monte Carlo (MC) simulations using alternating 
local Metropolis moves of single particles and collective cluster moves, specifically swaps 
of clusters between the configuration and its inversion at a randomly chosen point 
[35, 36]. Second, a 2D and a 3D event-driven molecular dynamics (MD) code solving 
Newton's equations for a set of particles interacting via a hard spheres potential [37, 38]. 
All simulations are performed in the NV ensemble with periodic boundary conditions; 
particle numbers are 16,000 in the MC simulations (4,000 for the larger densities 77 > 0.6 

\ If the random variables X, Y are independent Gamma variates with the same length scale of the 
Gamma distribution, then X/(X + Y) is a Beta variate. In the case of (3^' b , X and Y correspond to the 
smallest eigenvalue A m i n and the increment A max — A m j n , which are approximately Gamma distributed, 
but not independent variables. 
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Figure 3. Average anisotropy indices f3 for equilibrium ensembles of hard disks (left) 
and spheres (right) vs. packing fraction r\. In the limit rj — 0, the hard spheres 
anisotropies are identical to those of the Poisson process (solid bullets). 

In 2D (left-hand side), changes of slope at the transition densities (magnified in 
inset) are observed. In 3D (right-hand side) the average f3 displays a discontinuity at 
the first-order phase transition. The values of fl^' are very close to those of /3q'° and 
have been omitted for the sake of clarity. The inset shows metastable states in or close 
to the coexistence region; symbols are Monte Carlo results, lines molecular dynamics 
simulations. 



in 3D), and up to 256,000 in the MD 3D simulations. The starting configuration is an 
ordered triangular or fee lattice. 

The pair correlation functions of the resulting configurations are found to be in 
agreement with the Percus-Yevick approximation in the fluid phase [39]. The mean 
square displacement of the particles between two sampled configurations in the fluid 
phase is as large as the expectation for completely independent configurations, indicating 
a sufficiently equilibrated system. Away from the phase transition densities, both in 2D 
and 3D, the pair correlation functions from MC and MD are in agreement. 

4-1. Averages of Anisotropy Indices 

In similar fashion as with the Poisson data, a periodic Voronoi tessellation is computed 
for each configuration of particles, and Minkowski tensors and anisotropy indices are 
assigned to each particle. Fig. 3 shows the averages of (3® ,b for hard disks and hard 
spheres ensembles as a function of density. It is seen that increased order corresponds 
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Figure 4. Standard deviations of the f3 distribution for equilibrium hard disks (left) 
and spheres (right). A linear decrease as r) — > 77 ma x as been divided out. All data 
points calculated from Monte Carlo simulations. In the hard disks system, a{(3) is 
continuous across the phase transition (magnified in inset). The hard spheres system 
show a discontinuity. The values of /Sj' are very close to those of /3 ' and have been 
omitted for the sake of clarity. Solid bullets on the plot frames represent the limit of 
the Poisson point process rj = 0. 



to more isotropic Voronoi cells in both hard disks and hard spheres. 

We find MD and MC results to be in close agreement with each other except for in 
the vicinity of the hard spheres coexistence region. There, naturally, metastable states 
persist due to finite system sizes and simulation times, and the state of the system 
depends heavily on its history. The metastable fluid branch can be traced even beyond 
the fluid/solid coexistence region by choosing random initial conditions [40] instead 
of the usual fee initial conditions (see fig. 3, inset in righthand side). Given enough 
simulation time, it would eventually form fee nuclei and crystallize. The metastable 
fluid branch is distinctly more anisotropic on a local scale than the solid phase, again 
confirming the rule of thumb that more ordered systems are more locally isotropic. 

In both MC and MD simulations, a discontinuity in the average anisotropy index 
/i(/3"' 6 ) is observed at the hard spheres phase transition, with the solid phase being more 
isotropic locally. At the same time, the standard deviation of the distribution of 
drops discontinuously across the transition; in the solid phase, it goes to zero almost 
linearly, o~((3®' b ) oc (77 — ^ max ), where ^ max is the maximum global packing fraction. 

For large packing fractions, the system approaches the fee state with isotropic cells 
only (/i = 1, a = 0). For small packing fractions, the Poisson limit is recovered, and fi, 
a attain their respective values for the Poisson point process (bullets on the plot frame 
in figs. 3, 4). Remarkably, tensors W^ and Wq'° yield very similar anisotropy indices. 
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No explanation for this data collapse has been identified. The shape of the distribution 
of f3 values, seen in fig. 2, remains qualitatively similar and unimodal from the Poisson 
limit f] = to the close-packed limit; its average value has a discontinuity at the phase 
transition. The development of the distribution function with varying packing fraction 
can be seen in fig. 9 in the appendix. 

The situation in the hard disks system is less clear, and the nature of the phase 
transition has been the subject of continuous debate; current opinion favors a two- 
step transition according to the Kosterlitz-Thouless-Halperin-Nelson- Young scenario as 
opposed to a very weak transition of first order (see [41] and refs. therein). 

Mean anisotropy indices show a change of slope at the transition densities (fig. 3, 
left-hand). In the transition region, we use systems of 64.000 disks in a hexagonal 
simulation cell with periodic boundary conditions§. Wo find no significant differences 
between systems of 16.000 and 64.000 disks, given enough relaxation time. The 
curves in fig. 3, left-hand side, and second moments /x(/3 2 ) behave linearly as prescribed 
by the lever rule, with coexistence densities of ^0.70 and 0.72. Binder et al. [41] report 
closer values of 0.706 and 0.718. Full equilibration of the hard disks system in the 
transition region, however, remains difficult even with more advanced algorithms than 
the ones used in the present study [42]. 

4-2. Correlations among Anisotropy Indices 

We find that the four (in 2D; six in 3D) different anisotropy indices display a qualitatively 
similar behavior. This shows that the anisotropy of the Voronoi cells is a generic 
feature and does not depend on the particular aspect of anisotropy being analyzed. 
Shapes can easily be constructed where the different Minkowski tensors yield widely 
different anisotropy indices. However, such shapes do not generally occur in a Voronoi 
tessellation. 

This motivates the use of a single anisotropy index instead of the full set of 
Minkowski tensors. To justify this, we calculate the correlation coefficients corr(X, Y) := 
cov(X, Y)/(a(X)a(Y)), cov(X,Y) := n((X - n(X))(Y - fi(Y))), of the anisotropy 
indices; most pairings of (5^ ,h are correlated with corr > 0.8 over the whole packing 
fraction range of the hard spheres ensemble. Only corr^' , ft®' 2 ) gets as low as 0.7 for 
the smallest packing fractions. The same is true for the hard disks ensemble, corr > 0.8 
for most pairings; only corr(/3Q'°, ft®' 2 ) gets as low as 0.65 for the most dilute systems. 
It is thus justified, for the anisotropy analysis of the hard spheres ensemble, to focus on 

2 

a single convenient anisotropy index, for example the translation invariant p{ . In the 
remainder of the article, we drop the indices u, a, b off the /3. 

§ A square box was also used. As expected, it leads to artefacts for large packing fractions (starting from 
about r\ rts 0.75), and = 1,ct(/3) = is never attained as cells are strained by the incommensurate 
simulation box. 



Local Anisotropy of Fluids 



12 



4-3. Comparison to Isoperimetric Ratio 

Measures of asphericity have been previously used, e. g. the isoperimetric ratio (also 
called shape factor [9, 10]), defined by ( := {Wi) 2 (Wq)~ 1 / 'u 2 in two dimensions, and 
( := (Wi) 3 (W )~ 2 /u 3 in three dimensions, with the volume oo D of the D-dimensional 
unit sphere. Starr et al. use an asphericity index based on closest neighbors distance in 
the Voronoi tessellation [11]. 

The distinction between asphericity, meaning deviations of the shape from a sphere, 
and anisotropy, meaning that the object has different spatial extension or surface area 
in different spatial directions, is subtle. Clearly, a sphere's isoperimetric ratio can be 
increased by applying an undulation to the bounding surface or by deforming it into a 
cube. Both leave the anisotropy unchanged. On the other hand, ( also changes when 
the sphere is deformed into an ellipsoid, becoming anisotropic. The anisotropy indices (3 
clearly distinguish between asphericity and anisotropy. Also, the Minkowski tensors not 
only provide a convenient anisotropy index, but also identify explicitly, by means of the 
eigenvectors, the preferred directions for each cell. Accordingly, we find that ( is less 
correlated with the /3 indices than the set of /3 among themselves; corr(£, ranges 
between —0.6 and —0.2 for all (3® ,b both in two and three dimensions, meaning that 
the Minkowski tensors probe a different aspect of the cell shapes than the isoperimetric 
ratio alone. 

It is noteworthy that the isoperimetric ratio follows a double-peak distribution in 
the hard disks system [9]. The distribution of the f3 indices, however, is unimodal, see 
fig. 2; it changes continuously across the solid/fluid phase transition. 

5. Lennard-Jones fluid 

Configurations of Lennard-Jones (LJ) fluids with 100.000 particles in periodic boundary 
conditions are generated using a parallel molecular dynamics code [43, 44]. The LJ fluid 
is a true thermal system, the phase diagram is two-dimensional, and includes a critical 
point. The precise location of the critical point depends on the numerical treatment 
of the LJ potential, in particular the cutoff radius and cutoff corrections [45]. Here, a 
cutoff of 2.5r is chosen. 

For the pointlike LJ particles, we consider, as previously in fig. 1, the local 
normalized (particle number) density = l/(pV^), where p is the global particle 
number density, and the volume := Wq(K^) is the volume of Voronoi cell 
Thus, = 1. Fig. 5 shows two-parameter histograms of the densities and 

anisotropy indices 0^ Q {K^>) of Voronoi cells in LJ configurations. Configurations from 
four points in the LJ phase diagram are shown, and the corresponding locations indicated 
in the phase diagram sketch (inset). In the high-temperature limit, above the critical 
point, we find a continuous change of the probablity distribution from a triangular shape 
(subfigure a) as is also found for the ideal gas (fig. 1) to a ellipsoidal shape (subfigure 
b) as is found in the hard spheres system (also seen in fig. 1). Below the critical point, 
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Figure 5. Histograms of the Voronoi cell anisotropy and density at four points in 
the Lennard- Jones phase diagram (sketched in inset): a) vapor state b) liquid state c) 
critical composition d) coexistence of liquid and vapor. The systems consist of 100.000 
particles each. The x axis shows a local normalized density given by the reciprocal cell 
volume V in units of the global density p, while the y axis shows the anisotropy index 
/3q'° . A more complete survey of the phase diagram is found in the appendix (fig. 10). 



coexisting liquid and vapor phases can be clearly distinguished as two components of the 
probability distribution (subfigure d) . These components grow closer as the temperature 
is increased, and eventually merge at the critical point (subfigure c). 

6. Robustness of the Minkowski Tensors 

A simple model of a thermal solid is defined by putting the germs on an ideal lattice, 
and allowing each germ to deviate by a small distance e*w from its ideal site, as in 
the Einstein solid. We generated a point configuration by drawing, independently for 
each lattice site, the random displacement vector e*W from a three-dimensional Gaussian 
distribution. The set of germs is then {g^ — I W + £ W}, where are the fee lattice 
sites. The noise amplitude is given by the root mean square displacement (RMSD) 
of the particles, RMSD = a/ jj,(\e^\ 2 ); since the expectation value of e W vanishes, 
RMSD = a(e*W). Increasing the noise amplitude eventually destroys lattice order und 
yields a Poisson point pattern. This model neglects correlations between displacements 
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Noise amplitude RMSD/a 

Figure 6. Left: Anisotropy index for Wi' (K^) and Edwards' configurational tensor 
C^ l \ applied to a fee Einstein solid. The noise amplitude is quantified by the root 
mean square displacement RMSD of the germs from their ideal lattice sites. The 
Minkowski tensor W®' 2 (K^) is computed from the Voronoi tessellation of the lattice, 
is computed with the Delaunay neighborhood and distance cutoff neighborhoods 
of the order of the noise amplitude. It is seen that, even for vanishingly small levels of 
noise, of a single lattice site is not an isotropic tensor, while W 1 ' (K^) is. Tensor 
isotropy is defined via the ratio of eigenvalues, as in eq. (6). The x axis and cutoffs 
are in units of the fee nearest neighbor distance a. 

Right: Illustration of point configuration with a degenerate Delaunay triangulation. 
A small perturbation (not drawn to scale) lifts the degeneracy and breaks the isotropic 
4-simplices into anisotropic triangles. 

of different particles that necessarily exist in a real solid, and does allow for overlap if 
the particles have finite radius. 

In fig. 6, as the noise amplitude is increased, /z(/3°' 2 ) decreases continuously from 
unity, similarly as it does when decreasing packing fraction rj from the maximum density 
state in the hard spheres model (fig. 3). In the limit of strong noise, the anisotropy index 
of the Poisson point pattern is recovered (/x(/3°' 2 ) — > 0.457). 

6.1. Robustness against Noise 

The Minkowski tensors of Voronoi cells vary continuously when germs as dislocated. 
The topology of the Voronoi tessellation may change; however, the shape of the Voronoi 
cells is a continuous function of the germ point locations. For example, it is possible 
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that a small dislocation of a germ creates additional facets in the Voronoi tessellation; 
however, since these new faces are very small, and their contribution to the tensor of 
the cell is weighted with their area, only a small change in the tensor is effected. 

This continuity property is not necessarily guaranteed for other, similarly defined 
shape measures. For example, Edwards et al. [12, 46] defined a configurational tensor 



with some suitable definition for the set NN of nearest neighbors. NN may be defined 
using the Delaunay triangulation, or a fixed cutoff radius, or physical contacts. If, 
under a small dislocation of a germ, the neighborhood NN changes, a comparatively 
large change of the tensor results, because the bond vector fW — r ^ to the lost or 
newly gained neighbor tends to be larger than the other bond vectors. Consequently, 
changes discontinuously|| . 

The discontinuous behavior of is especially pronounced for near-degenerate 
configurations, for example the fee lattice with vanishingly small noise amplitude. As 
fig. 6 shows, is not an isotropic tensor at zero noise, even though a fee lattice site has 
12 symmetrically arranged nearest neighbors. This is due to numerical roundoff errors. 
For simplicity, we discuss the two-dimensional case, which exhibits the same conceptual 
problem. Consider the perfect square lattice in fig. 6. The Delaunay simplices are not, 
as in a general point pattern, triangles, but degenerate to squares. computed from 
such simplices is an isotropic tensor. However, tiny amounts of numerical inaccuracy or 
noise lift the degeneracy for some of the squares and cause the squares to break up into 
triangles [29]. computed from the triangles is no longer an isotropic tensor. Due to 
the same effect in three dimensions, with the Delaunay neighborhood is not a good 
measure of anisotropy. 

2 

Note that the Minkowski tensor W 1 ' is similar to the Edwards tensor as the normal 
directions of Voronoi facets are precisely the directions of the bond vectors via the 
Voronoi-Delaunay duality. However, the weighting with the facet area in eq. (5) ensures 
that the tensor is a continuous function of the germ point locations. 

6.2. Discrimination of Einstein Solid and Hard Spheres Equilibrium Solid 

As shown above, the mean anisotropy index //(/?) of the Einstein solid decreases 
with increasing noise amplitude, similar to the hard spheres solid. However, the two 
point patterns are clearly structurally different, as the hard spheres system includes 
correlations between individual particles. 

To see this, it is again useful to study the correlation between local normalized 
density d and local anisotropy (3. Figure 7 shows that is negatively correlated with 
j3(K^) for the Einstein solid (blue dashed line), while the two variables are positively 

|| Note that Edwards et al. use to describe physical contact points in granular matter, and for 
that purpose, must change discontinuously. However, as a consequence, is not suitable to 
characterize local anisotropy in a robust way. 
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Figure 7. Comparison of average anisotropy index /3°' 2 as a function of local 
normalized densities d. The Voronoi cells are binned in density bins, and an average 
of their anisotropy indices is computed separately for each bin. The error bars are 
estimated from standard deviations in each bin. The plot shows that the correlation 
of anisotropy index and cell volume distinguishes between the Einstein solid and the 
equilibrium hard spheres structure (in the solid phase), even though parameters were 
chose such that the global averages fi{fiy ) are equal. 



correlated for the hard spheres solid (red solid line). Positive correlation is the rule for 
interacting systems such as hard spheres, and has also been found for jammed bead 
packs [4], while a negative correlation has been found for foam models. 

This can be understood considering that, for the hard spheres solid, a higher density 
implies that the cell shape is closer to the shape of the particle since particles cannot 
overlap. On the other hand, for the Einstein solid, no minimum distance is enforced, 
and high densities are the result of one germ intruding into the domain of another, 
creating two oblate Voronoi cells. 

Conclusions 

This study shows how Minkowski tensors can be used to characterize point patterns, 
and in particular their local anisotropy. A simple, linear algorithm to compute the shape 
tensors of a polyhedron was given elsewhere [25, 26]. 

Particle configurations from a number of fluid models are analyzed using the 
Minkowski tensor framework, and reference values for these idealized systems are 
given. In each of the systems, the Minkowski tensor analysis is able to quantify the 
morphological properties of the system in a robust way. The continuity properties 
guarantee that small shifts in the germ positions do not influence the result in an undue 

2 

manner, as is the case for a naively defined "bond tensor" . In particular, W 1 ' is found 
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to include the correct weighting prefactors of the bond vectors to ensure continuity. 

The anisotropy analysis is found to be sensitive to the phase behavior of the fluids, 
and the average anisotropy index may be used to distinguish collective states of the fluid. 
More research is needed in order to link order parameters based on Minkowski tensors 
with more established order parameters like the bond-angle [5] order parameter. For the 
hard spheres system, exact relations between the equations of state and the geometry 
of the available space, in particular its surface area and volume, are known to exist 
[47, 48, 49, 50]. Similar relations are not yet known for Minkowski tensors of available 
space. However, a density functional theory based on Minkowski tensors of convex 
particles has been developed, and provides accurate approximations for the equations of 
state [51, 52, 53]. Finally, the application of the Minkowski tensors to a comprehensive 
set of experimental data is desirable. For the latter purpose, the Minkowski tensor 
software may be freely downloaded from the institute homepage at 

http: //www.theoriel .physik.uni-erlangen.de/ 
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Figure 8. Distribution of scalar cell properties in a Poisson-Voronoi tessellation: 
The probability density is logarithmized, and its argument replaced by the reduced 
quantity w v :— InWt, — ^,{h\W u ). In this plot, a log-normal distribution would 
appear as a parabola with the apex on the w — line. The data is generated from 
5 x 10 7 (2D) and 10 7 (3D) random Voronoi cells. Generalized Gamma distributions 
oc (W,,)^" -1 ) exp(— {W u /9) q ) give good fits with the parameters specified in the plot 
key. {6 sets the length scale and is fixed by the normalization condition fi(w u ) = 0.) 
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Figure 9. Histograms of anisotropy index /3 ' and packing fraction rj of Voronoi 
cells in the hard spheres ensemble. The x axis of each box shows the packing 
fraction jjW g [0;0.75] of individual cells, the y axis shows the anisotropy parameter 
/3 ' [K^>) s [0;1]. Darker gray corresponds to larger probablity. The distribution 
changes characteristically from a triangular shape (77 < 0.2) to a ellipsoidal distribution 
for higher densities. In the limit 77 -4- n/y/TE « 0.74, the distribution becomes a Dirac 
distribution f(r],P) = 5{q - ir/Vts)6(l - /?). 
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Figure 10. Histograms of /3q'° and d for Lennard- Jones configurations. The x axis 
of each box shows local particle number density dW g [0;4.5] of individual Voronoi 
cells, the y axis shows the anisotropy parameter /3 ' (K^) € [0; 1]. The critical point 
would correspond to the configuration at T = 1.1 and p = 0.26; below, vapor/liquid 
coexistence is found; above, a single fluid phase exists and the probability distribution 
only has a single component. 



